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Abstract 

We construct quantum circuits for solving one-dimensional Schrodinger equations. Simulations 
of three typical examples, i.e., harmonic oscillator, square- well and Coulomb potential, show that 
reasonable results can be obtained with eight qubits. Our simulations show that simple quantum 
circuits can solve the standard quantum mechanical problems. 
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I. INTRODUCTION 



Quantum computers have been one of the most growing fields in computational physics 
last two decades. Since Feynman suggested that a quantum computer could possibly simu- 
late quantum systems more efficiently than a classical one a large amount of work has 
been devoted to quantum algorithms and their experimental realizations. This is because 
the quantum register (qubits) can store data in the superposition of quantum states and 
operations of them can be executed in parallel, which results in an exponential reduction of 
the computation time and memories. Among powerful applications of quantum algorithms 
are Shor's factoring integers Q and Grover's searching databases {^J. 

Since few quantum circuits are universal, i.e., any unitary operation on qubits can be 
constructed by those universal quantum gates, it is, in principle, possible to make appro- 
priate quantum circuits for calculating classical functions 0, 0] . However, the efficiency of 
simulations with quantum circuits is very much dependent on the dynamical system under 
consideration. Therefore we must find an efficient way of describing the system and an 
efficient quantum simulation algorithm. 

So far, there have been proposed several quantum algorithms for simulating quantum 
mechanical systems. Simulations of many body system have been reported, i.e., lattice-gas 
[(J, Heisenberg model 0, 0], pairing model jj 10] and Hubbard model [lO, [ill, Q. Also 
quantum computations are expected to provide polynomial-time simulation of chemical dy- 



namics 



13, 



14j . These systems are suitable for the simulation with quantum algorithm, since 



their quantum states are naturally represented by qubits, i.e., |0)/|1) of a qubit corresponds 
to the eigenstate of the number operator in the second quantized formalism, or the up/down 
state of a spin, for example. On the other hand, there have been few quantum simulations 
for particles in real space, although a general algorithm was developed by Zalka 15] and 
Wiesner [la]. Recently, Benenti and Strini 17 1 simulated time-evolution of a Gaussian wave 
packet in the harmonic oscillator potential, and Oh [l_8| calculated the ground state energy 
of a displaced harmonic oscillator and a quartic anharmonic oscillator. 

The purpose of this paper is to provide concrete examples of explicit simulation of the one- 
dimensional Schrodinger equations of typical potentials, i.e., harmonic oscillator, square- well 
and Coulomb potential. We will explicitly construct the quantum circuits for the calculation 
of the eigenvalues and eigenstates of those Schrodinger equations. In section 2, we will 
describe how to make quantum circuits for these three examples. Several simulations will be 
reported in section 3, in which outputs are compared with exact values. Section 4 is devoted 
to a summary. 
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II. QUANTUM ALGORITHM FOR SOLVING SCHRODINGER EQUATIONS 



In 

tions 



;his section, we will briefly review the quantum algorithm to solve Schrodinger equa- 

mQQ. 



IS). 



A. Time-evolution of the quantum state 

Let us consider the case where a single particle is moving on a line under the potential 
V(x). The one-dimensional Schrodinger equation is 



H\ 



2m 



£- + v{x) \^) = ih—\^) . (l) 



dt 



Hereafter, we set the mass m = 1 and Plank's constant h = 1 for simplicity. In the case of 
time-independent Hamiltonian, the formal solution of Eq.flT]) is 



\m) = u(t)\m) = e-nm) , w 

where U(t) = e~ lHt is the unitary operator of time-evolution. In order to calculate the 
time-evolution, firstly time-interval t is divided into n steps, i.e., t = nAt, and then each 
step is approximated by the second-order Trotter formula as 

e -iHAt = e -i(K+V)At = e -iV At/2 e -iK At e -iV At/2 + ^ ^ 

where K = p 2 /2 is the kinetic operator. While it is straightforward to calculate e~ lVAt ^ 2 in 
the coordinate basis \x), it is preferable to calculate e~ lKAt in the momentum basis \p). The 
transformation of the basis is defined by 

/oo ^oo 
dx\x)(x\p) = dx e 2wipx \x) , (4a) 
-co J — OO 

/OO POO 
dp\p)(p\x)= dp e~ 2mpx \p) , (4b) 
-oo J — oo 

where we have used the convention for the later convenience. The coordinate representation 
of the state vector lib) is 



dx\x)(x\ip) = / dxijj(x)\x) , ip(x) = (x\i/j) , (5) 

) J — oo 

and, the Fourier transformation of the wave function is 

/oo 
dxe~ 27ripx ij(x) = Ul T t/j(x) , (6a) 
-oo 

/"OO 

ip(x) = (x\if>) = / dpe 2nipx ijj(p) = U FT i){p) . (6b) 
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Therefore, x-representation of the wave function \ip(t)) = U(t)\ip) is 



/ x \ e -iVAt/2 e -iKAt e -iVAt/2i 



oo 



= e -m*)M/2 / dpe 2 mpXe -iK( P )M / dx ' e -^' e -iV( X ')At/2^ 
J —oo J —oo 

= e" iy(:c)A * /2 f/ FT e- 4K(p)Ai 4 T e- iV ' (a; ' )A * /2 ^(x / ) . (7) 

B. Discretization of the coordinate and quantum Fourier transformation 

We are interested in the bound state where the wave function ip(x) is localized in some 
finite region. The wave function ip(x) can be approximated on appropriate mesh points {xk} 
in this region as 

|V) = 5>(zk)kk> • (8) 
k 

In order to carry out the Fourier transformation Eq.plh we will employ quantum Fourier 
transformation (QFT). 

The QFT is the unitary operation which transforms the basis {|0), . . . \N— 1}} to the 
new basis {|0), |1), . . . , \N— 1)} such that 

N-l 

Uqft : \j) - |j> = ^=E <™>" N \k) , (9a) 



N k =o 



Uqft ■ \k) - l*> = "4 Yl e- 2 ^ N \j) . (9b) 



3=0 



By comparing with Eqs.(j4]), the bases \k) and \j) are identified with the coordinate basis 
\xk) and the momentum basis \pj) respectively. Then inverse QFT of Eq.(jBJ) is 

iv>> = E^)i fc > = Y,^4^ e ~ 2nl3k/N & = • (io) 



Therefore the inverse QFT changes the x-representation of ip(x) to the p-representation ip(p) 



15 



16|]. 



By making suitably scaling and shifting the coordinate, we will choose the x-space interval 
[—1/2, 1/2] for simplicity, and iV equally spaced mesh points, i.e., x^ = k/N — 1/2, (k = 
0,1, . . . , N — 1). Accordingly, the p-space mesh points are taken as pj = 2n(j — N/2), (j = 
0,1,..., N-l). In this case, QFT Eqs.® are 



\j) = _Lj2e 2mU - N/m/N ~ 1/2) \k) , (11a) 

fc=0 

1 N ^ 

\k) = ^Y, e ~ 2 ^~ NmklN ~ m \~3) • (Hb) 



N ■ n 

3=0 
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The phase factor becomes e 2m i k l N ' e m ( k +i) e 2mN / A . The constant factor e 2mN ^ 4 can be ab- 
sorbed in the bases. Redefining the new bases 

\ k >) = e~ mk \k) = (-l) k \k), = e^\j) = (-iy\j) , (12) 

the standard QFT Eqs.([9]) can be satisfied. Therefore, in executing the practical calculation 
with the standard QFT, the wave function should also be redefined as 

|V) = X>(zfc)l*> = E^ Xfc )l fc '> ' = W^k) • ( 13 ) 

k k 

The distribution of mesh points {xk} described above is not exactly symmetric with 
respect to x = 0. This may cause some numerical inconvenience for symmetric potentials. 
The mesh point Xk = is also not suitable for Coulomb potential. Therefore it is convenient 
to prepare another distribution which is exactly symmetric and does not contain the point 
x — 0. They are 

Xfc = ^-G"2]v)' ^ = 27r [ J -(f-^)] • (M) 
In this case, the Eqs.( !TT]) become 

N-l 

N 



\j) = 1 ^ e 2-(i-iV/2+l/2)( fe /iV-l/2+l/2iV)| A . ) ^ (15a) 

v N k=0 
N-l 

= 1 y e -2iri(j-N/2+l/2)(k/N-l/2+l/2N) jj^ _ 



^^0 



Then, with new bases 

|^ = ^(1/2-1/2^^ ^ = e 2^(l/2_l/2iV) i |J ) ) (16) 

the standard QFT is satisfied. The wave function is accordingly redefined as 

${z h ) = e™^-y™^{x k ) . (17) 
This distribution will be employed in the next section. 



C. Phase estimation and eigenfunction 



If we take the initial state | -0(0) ) to be an eigenstate \uk) of the Hamiltonian H with 
eigenvalue E k , i.e., H\uk) = Ek\uk), then U(t)\u^) = e~ lEkt \uk) and the energy Ek can be 
calculated by the phase estimation algorithm 20j. Since the phase estimation algorithm 
finds the eigenvalue e 2m ^ of the unitary operator U(t), the energy eigenvalue Ek is given by 



Ek 



-2n(j)/t, (0 < <j) < 1) 



(18) 
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In order to make the energy eigenvalue negative, it is necessary to shift the Hamiltonian by 
an appropriate constant value. One should also choose the evolution time t such that the 
searched energy range is (E max — E min ) = 2-n/t. Since the phase estimation algorithm gives 
us the same E^ periodically, we should be careful about the situation where different energy 
states may contribute the same energy phase. 

The phase estimation algorithm consists of two kinds of registers, i.e., the first register is 
work qubits for storing the phase of the unitary operator U(t), and the second register is the 
simulation qubits for representing the quantum state. The total state is a tensor product of 
work qubits and simulation qubits. In general case, the initial state \ip(0)) is written by the 
superposition of eigenstates of H as 



Thus, the total state is also the superposition of the tensor products, and one can find the 
coefficient by the projection operator of the corresponding work qubits 0, [l^]. In order 
to execute efficient simulations, the initial state should be prepared in such a way that it 
has an appreciable overlap with the eigenstate |tt&) which we are searching. 

III. SIMULATIONS OF TYPICAL EXAMPLES 

In this section, we will show the practical way of constructing quantum circuits for three 
typical examples, i.e., harmonic oscillator potential, square- well potential, and Coulomb 
potential (S'-wave). The range of the coordinate x is fixed to [—1/2, 1/2], and we will choose 
the strength of the potential such that the wave function is localized in this range. We 
will employ w work qubits (first register) and s simulation qubits (second register). Thus 
the dimensions of the work and simulation bases are N w = 2 W and N s = 2 s respectively. 
Total number of qubits is q = w + s and the dimension is N q = 2 q = N W N S . Then, we 
will prepare equally spaced N s mesh points for —1/2 < x < 1/2, and N w energy points for 
2ii/t = (E max - E min ) with energy step size AE = 2ix/t/N w . 

In the practical calculations of the following examples, we set w = s = 4. The typical 
energy scale is 10 2 , and simulations give good convergence with divided time interval At = 
t/n c± 10~ 3 . We have carried out several calculations with more qubits and time-steps, 
and certainly obtained improved results, although the qualitative features remain the same. 
Therefore we will show the results of simulations with w = s = 4, which can be executed 
within reasonable computer resources. In our experience of numerical calculations on the 
ordinary (classical) computer, 2 4 = 16 mesh points or bases are sufficient to obtain ground 
and a few excited states in one- dimensional potential. So it is expected that simulations 
with w = s = 4 could give us outputs with more or less similar accuracy. 

The quantum circuits for the quantum Fourier transformation (QFT) and the phase 




(19) 



k 
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estimation are well known and detailed descriptions are given in Ref. 19| for example. Thus 



we will not repeat the explanation of these circuits. In Ref. 19|, one can also find how 
efficient is the quantum simulation algorithm. 

In the following subsections, we will explicitly construct quantum circuits of the time- 
evolution operator U(t) = e~ tHAt , execute simulations with appropriate initial states, and 
compare the outputs with exact values. For these examples, quantum circuits can be con- 
structed only by single- and two-qubit operators. Furthermore, ancillary qubits calculating 
the potential term are not necessary. The phase-evolution due to the potential term is 
implemented directly in the time-evolution circuit. 

A. Kinetic energy term 

Let us begin with the quantum circuit of the common kinetic energy term e~ iKAt . The 
time-evolution operator of the kinetic term is 

e- iKAt \ Pj ) = e-'^lpj) = e ia{ ^- l 2?\ pj ) , a = -(2nN s ) 2 At/2 . (20) 

The integer j is represented by the binary form as, 

s 

j = J2jn2 s ~ n = ji? ; + j->T 2 + . . . + j s 2° = j lj2 . . . j 8 (binary) . (21) 

71=1 

Thus j/N s G [0, 1] is the binary fraction 

s 

j/ N * = J2^ 2 ' n = J >' 2 ' • / -- 2 + • • • • -I- ' = Q W2 ■ ■ -j s (binary) . (22) 

n=l 

Therefore 

(j/N s - 1/2) 2 = (J> n 2- n - 1/2) 2 = (J2 3n2- n f - + 1/4 . (23) 

n n n 

The computational basis \pj) is the direct product of s qubits 

\Pj) = \j1j2 ■ ■ ■ Js) = \ji) ® |j 2 ) . . . ® \js) ■ (24) 

The last term of Eq. (|23p simply multiplies a constant factor e %a / A on one of the simulation 
qubits, \jx) for example. The operation of the second term is 

e — E„i-2-| jlj2 . . _ jg) = (S?) e -^ 2 -"|j„) =<S$R{-aQ- n )\j n ) , (25) 



where R(9) is the single-qubit operator rotating the phase of |1) by 6, i.e., 

rm = I 1 I . (20) 



e 



w 



The first term of Eq. (!23l) is written by 



<(E m j m 2- m )(e„j~2-") 



exp(za ^ j^- 2 - + 2za £ J m J„2- m - n ) 



(27) 



rn=£n 



Since j 2 = j n) e 40 - 7 ™ 2 2 " is a single-qubit operator given by R(a2 2n ). On the other hand, 



Jmjr 



jm = Or j n = 

1 Jm jn 1 i 



(28) 



2iaj m j n 2 



is given by a two-qubit operator A as 



.4 



/ 1 





\ 


1 







1 







\0 


e 2ia2- m - 





(29) 



This two-qubit operator acting on \j m j n ) can be represented by the controlled-?/ {CU) 
operation shown in Fig.l. 



\jr. 



\3n) 











u 





FIG 1: Quantum circuit CU 



In the CU circuit, the single-qubit operator U is applied to the target qubit when the con- 
trol qubit is set to |1). In this case, U = R(2a2~ m ~ n ). For s simulation qubits, the number 
of two-qubit operator is S C*2 = s(s + l)/2, and one should apply CU gates successively. 



B. Harmonic oscillator potential 

The Hamiltonian of the harmonic oscillator is 

H = \p 2 + y* 2 . (30) 

Since the potential term is the same quadratic form as the kinetic term, the quantum circuit 
of the time-evolution e~* yA '/ 2 is the same as the kinetic energy term, except that the strength 
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a = — (2nN s ) 2 At/2 is replaced by (3 = ~uj 2 At/A. Then, the calculation of the time-evolution 
operator is 

e -iVAt/2 e -iKAt e -iVAt/2^ = ^WAt^^-iKAt^^-iVAt^ ^ 

We have chosen the potential strength parameter uj such that both x- and p-space wave 
functions ip(x) and ip(p) are well localized in the chosen finite interval and transformed 
accurately by the QFT. In practice, optimal parameter uj is given by 

uj/2 = (2nN s ) 2 /(2uj), uj = 2nN s . (32) 

For s = 4, uj ~ 100.53 and we fixed uj = 100 in the following calculations. 

We will show the probability spectrum \ck\ 2 of Eq. (1191) as a function of the energy E in 
Figs. 2 for the initial states ipo(x) = ^-^l 2 and ipo{x) = xe~^ x ' 2//2 (the normalization factor 
will be omitted hereafter). Since these initial states are exact eigenstates, the outputs are 
good check for the simulation. The energy spectrum of Fig. 2 (a) clearly shows the sharp peak 
around the exact value E = uj/2 = 50. The numerical value is |c| 2 = 0.915 at E = 52.4. 
Fig. 2(b) shows the result of the first excited state and output of simulation is |c| 2 = 0.699 
at E — 157. These examples show that our simulations work fairly well. 




20 40 60 80 100 120 25 50 75 100 125 150 175 200 



(a) (b) 



FIG. 2: Probability spectrum, (a) if>o(x) = e"^' 2 / 2 . Parameters t = 0.045, n = 30. (b) 
ipo{x) = xe'^ 2 / 2 . Parameters t = 0.03, n = 20. 

Fig.3 shows the result of the initial state tpo(%) — x 2 e~ U)x2 / 2 . 

Although this is not the exact eigenstate, it is a superposition of the ground state <f>o(x) 
and the second excited state fafe), be., 

^o(x) oc y|</>o(z) + ]j^M x ) ■ ( 33 ) 

Therefore the energy spectrum shows two peaks around E ~ 50, 250, and the ratio of the 
probability is roughly 1 : 2, as is expected. 
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0.6 
0.4 
0.2 



50 100 150 200 250 300 



FIG. 3: Probability spectrum with ipo(x) = x 2 e ^ I 2 . Parameters t = 0.02, n = 20. 

Figs. 4 show the result of the initial state ipo(x) = 1/ cosh 2 (20:r). This state is a superpo- 
sition of even eigenstates. Fig.4(a) shows the peak at E — 52.4 with probability |c| 2 = 0.61, 
which corresponds to the ground state component. The exact overlap value is |c| 2 = 0.68, 
which is in good agreement. The small bump around E ~ 110 may come from higher ex- 
cited states. Fig. 4(b) shows the spectrum where the searched energy range is extended up to 
E = 300, and the second peak at E ~ 250 can be seen clearly . The ratio of the probability 
is about 4:1, which is also in good agreement with the exact value. 




20 40 60 80 100 120 50 100 150 200 250 300 



(a) (b) 

FIG. 4: Probability spectrum with ipo(x) = 1/ cosh 2 (20x). (a) Parameters t = 0.045, n = 30. (b) 
Parameters t = 0.02, n = 20. 

Figs. 5 show the projected eigenfunction corresponding to the peak at energy E = 52.4 of 
Fig.2(a). 

The wave function is normalized to be real at x — 0. The solid line is the exact eigen- 
function (f>o(x). Since the initial state is an exact eigenstate, the good agreement means that 
the QFT and the phase estimation algorithm work properly. The magnitude of imaginary 
part shows the inaccuracy of this simulation. 

One may wonder the outcome if the initial state is chosen randomly, which might cor- 
respond to ab initio calculation. Fig. 6 shows the average of the outputs of 10 random 
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(a) 



(b) 



FIG. 5: Projected eigenfunction of the ground state. Solid line shows the exact wave function, 
(a) Real part, (b) Imaginary part. 

initial states. There are three broad peaks corresponding to the exact energy values. This 
simulation shows that the initial state should be prepared carefully. 

0.25 

0.2 
.15 

. 1 
.05 

50 100 150 200 250 300 




FIG. 6: Average of 10 random initial states. 



The example of the harmonic oscillator shows that the quadratic potential can be con- 
structed by single- and two-qubit operators. One can readily understand that, for general 
n-th order polynomial potential, the quantum circuits are given by 2-, 3-,. . . ,n-qubit oper- 
ators, i.e., CU, CCU, . . . , C n ~ l U gates. 



C. Square- well potential 



The Hamiltonian of the square-well potential is 

H = \p 2 + V(x) , (34) 
where the potential energy V(x) is given by 

( -V„ 

V{x) = { ' ' (35) 



f -V 


\x\ 




\x\ 
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We will fix the potential strength V = 100 and the range a = 1/4 hereafter. Since the 
mesh points are distributed in [—1/2, 1/2], we choose these parameters such that the wave 
function is localized in this region. And it also makes the quantum circuit very simple, 
although the modification for general case is straightforward. 

Using the binary representation of k — j±j 2 ■ ■ -j s , the coordinate x k is written by 



k/N 8 - 1/2 = ^2j n 2~ n - 1/2. 



(36) 



71=1 



Thus, the first two qubits determine the position of x, i.e., 

-l/2<x< -1/4 for 3l = 0, j 2 = , 

-1/4 < x < for 3l = 0, j 2 = 1 , 

< x < 1/4 for ji = 1, j 2 = , 

1/4 < x < 1/2 for ^ = 1, j 2 = 1 . 

Therefore the potential energy becomes 



(37a) 
(37b) 
(37c) 
(37d) 



V 



ji = J2 = or j\ =32 = 1 
-V ji = 0, j 2 = 1 or j\ = 1, j 2 = 



(38) 



The time-evolution of the potential term e lVAt / 2 can be expressed by the two-qubit operator 
working only on the first two qubits \jij 2 ) as 



B = 



f I o o o\ 

e iy ° A */ 2 
e iy ° A */ 2 
\0 1 J 

This circuit B is constructed with a single-qubit operator 



(39) 



B 



e iV At/2 q 
q e iV At/2 



(40) 



and X-operator (NOT-circuit) which exchanges the coefficients of a single-qubit as 



X = 




(41) 



and controlled-?/ operator. Fig. 7 shows the quantum circuit executing two-qubit operator 
B, where the empty circle indicates that the operation is applied on the target qubit when 
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X 



X 



B 



X 



X 



-e — 6 — ©- 



-e- 



x 



x 



-e- 



FIG. 7: Quantum circuit B 




FIG. 8: Probability spectrum, (a) ^{x) = e 10x . Parameters t = 0.06, n = 50. (b) ipo(x) = 
xe ~Wx 2 _ p aram eters i = 0.06, n = 50. 

the control qubit is set to |0). The symbol © shows the X-operator (NOT-circuit). 

Figs. 8 show the probability spectrum |c fe | 2 as a function of the energy E for initial states 
ipo(x) = e~ Wx2 (even state) and tpo(x) = xe~ 10x2 (odd state) respectively. 

The exact energy levels are E = -88.12, E x = -54.05 and E 2 = -7.005. Fig.8(a) shows 
a sharp peak at E ~ —85 corresponding to the ground state, while a small bump at E ~ — 7 
corresponds to the second excited state. Fig. 8(b) also shows a sharp peak at E ~ —55 
corresponding to the first excited state. 

Figs. 9 shows the projected wave function corresponding to the energy E = —85.08 of 
Fig.8(a). 

The phase of the wave function is normalized as lm(^(0)) = 0. The exact wave function is 
shown by a solid line. The agreement is not so good as compared with the harmonic oscillator 
case. The wave function is slightly asymmetric, i.e., shifted to the negative direction, and 
also shows a strange behavior at \x\ ~ 0.3. The mixture of the imaginary part is not small, 
which clearly shows that the simulation has some problems. This is mainly caused by the 
fact that the potential is not exactly symmetric. This is due to the asymmetric distribution 
of the mesh points {xk}. Namely, at the boundaries of the potential- well, the strength is 
V = —Vq at x — —1/4 corresponding to \j1j2) = 1 1 ) , while V = at x = 1/4 corresponding 
to \j1j2) = 1 11)- Therefore the potential- well is negatively shifted by 5x = 1/2 5 in this case. 
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FIG. 9: Projected eigenfunction of the ground state. Solid line shows the exact wave function, 
(a) Real part, (b) Imaginary part. 

Another reason may be due to the sharp change of the potential at the boundary. 

In order to improve the simulation, we have employed the symmetric distribution of the 
mesh points given by Eq. fll4p . which also makes the potential exactly symmetric. The result 
is shown in Fig. 10 and Figs. 11. 

0.8r 

0.7 - 
0.6- 
0.5- 
0.4 - 
0.3 - 
0.2- 
0.1 - 

-100 -80 -60 -40 -20 




FIG. 10: Probability spectrum with iI)q(x) = e 10x . Parameters t = 0.06, n = 50. Mesh points 
are symmetrically distributed. 

Now the probability shows the more pronounced peak at E = —85.08. The phase of the 
wave function is set to real at \x\ = 1/2 5 . Figs. 11 show that the agreement of the calculated 
wave function with the exact values becomes much better. 

Fig. 12 shows the average result of 10 random initial states. Although two lowest states 
(E ~ —88, —54) may be seen as broad peaks, the third state (E ~ —7) cannot be resolved, 
and small fractions of many eigenstates seem fill over wide energy range. 
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FIG. 11: Projected eigenfunction of the ground state with symmetric mesh points. Solid line 
shows the exact wave function, (a) Real part, (b) Imaginary part. 

0.2 - 
0.15- 




FIG. 12: Average of 10 random initial states. 

D. Coulomb potential 

The Hamiltonian of the Coulomb potential is 

H= 1 -p"~-, («>0). (42) 
2 r 

The Schrddinger equation is reduced to one-dimensional equation in the case of S"-wave. Thus 
the solution ip(r) is given by ip(r) = rR (r) (r > 0), where Ro{r) is the radial part of the 
S'-wave Coulomb wave function. If the potential V(x) = —k/\x\ is defined in — oo < x < oo, 
the energy eigenvalues are doubly degenerate except for the ground state. We will take the 
odd wave function ip(x) = xRq(\x\) by setting odd initial states, since it is smooth at x = 0. 
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The problems of the one-dimensional Coulomb potential have been discussed in Refs. 
in detail. 

The construction of the quantum circuit of the Coulomb potential is not straightforward, 
since it is necessary to express the inverse of the binary fraction. We have made a simple 
expression in the following way. 
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Let < x < 1 be expressed as the binary fraction as 

N 



X = 

k=l 



Jk? k = O.J1J2 ■■■3n (binary) . (43) 



We will find the formula y = 1/x in terms of ji,j2, ■ ■ ■ ,3n- The first bit j\ can be set 
ji — 1 (1/2 < x < 1). In the case of j± = (0 < x < 1/2), one can shift the binary 
expression by an appropriate power of 2. Since y = 1/x > 1, y can be expressed by the 
power series of 1/2 as 

y = l + J2 a i 2 ~ e - (44) 

£>1 

Note that the integer coefficients are not necessarily or 1. In fact we find that a/s are 
small integers and the Eq. (14*41) converges rapidly. 
The coefficients are determined by the equation 



N 

x x y 

k=l 1 ! 

iv 



(E^ fe )( 1 + E a ^ 2 " 



k=l k,£ 

= 1 . (45) 

Since 1 = l/2 + l/2 2 + l/2 3 + ... = 0.111... in the binary fraction, one can obtain the 
equations which determine the coefficients ag recursively, 

3m + h a t = 1, m = 2, 3, . . . . (46) 

k+l=m 

Up to = 7, a/s are expressed as follows, 

(47a) 
(47b) 
(47c) 
(47d) 
(47e) 

+2j 2 j 6 + 2j 3 j4 + 2j 3 j5 - 6J2J3J4 • (47f) 

In the case ji = 0, one can obtain similar expressions by shifting j m — > j m+ i and multiplying 
by 2. 

The Coulomb potential is an even function and it has a singular point x = 0. Therefore 
the exactly symmetric mesh points of Eq. ffMj) is suitable. For the case of simulation qubits 



ai = 


1 


-32 , 






a 2 = 


1 


-33 , 






a 3 = 


1 


-h - 


33 


- 34 + 2j 2 j3 , 


04 = 


1 


-34 - 


35 


- 3233 + 2j 2 j4 , 


a 5 = 


1 


-32 - 


34 


- js ~ 3d - 3234 + 2j 2 j5 + 2j 3 34 


a 6 = 


1 


-33 - 


35 


- 36 - 37 + 3233 + 3j2j4 - 3235 
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s = 4, mesh points are explicitly given by 



x k = x - (1/2 - 1/2 5 ) = O.JU2J3J4 - 0.1 + 0.00001 (binary) . (48) 
The potential is proportional to the inverse of the absolute value \xk\, which is given by 

{O.OJ2J3J4I (binary) for j x = 1 
(49) 
O-Oj^l 1 (binary) for ji = , 

where = 1 — j m is the bit-flip of j m . Note that we can formally set j'5 = 1 for both cases. 
Thus, for Xk < 0, one should apply bit-flip operation before executing the time-evolution 
operator. 

The time-evolution operator of the potential term e~ lVAt ^ 2 can be constructed recursively 
depending on whether the qubit is |0) or |1). Defining the projection operator Pq (Pi) to 
the qubit |0) (|1)), 

' 1 o\ f o\ 

, Pi = , (50) 

J \o 1 / 

the matrix U\ (2 3 x 2 3 ) corresponding to x k > is given by 

Ui = P0U2) ®U 2 + P 1 (j 2 ) ® e ^^0-3.i4)At/2 ) (51a) 

^2 = PoO's) ® ^3 + PiUs) ® e"^^)^/ 2 , (51b) 

with 



^ 2 (j3, J 4 ) = 2{1 + (1 - j 3 )2- 1 + (1 - J 4 )2- 2 + (-J3 - J4 + 2/:,/;)2 :! 

+ (2j 3 - J3J 4 )2" 4 + (-2j 3 + 2j 4 )2- 5 + (1 + 3j 3 + J4 - 5j3j4)2- 6 } , (52a) 
V 3 (j 4 ) = 2 2 {1 + (1 - j 4 )2- 1 + j 4 2- 3 + (1 - j 4 )2- 4 + (1 - j 4 )2" 5 + j 4 2- 6 } , (52b) 
V 4 = 2 3 (1 + 2~ 2 + 2~ 4 + 2" 6 ) , (52c) 
C/ 4 = 2 4 (1 + 2- 1 + 2- 2 + 2- 3 + 2- 4 + 2~ 5 + 2- 6 ) . (52d) 

These formulas can be obtained by appropriately modifying the basic formula Eqs. fT4T|) . The 
time-evolution operators are single- or two-qubit operators, and can be constructed in the 
same way as the harmonic oscillator case. 

The simulations are carried out with a strength parameter k = 10. The accuracy of our 
approximation of the Coulomb potential with s = 4 simulation qubits is within 1.6%, which 
might be sufficient for simulations. 

Fig. 13 shows the probability spectrum as a function of the energy E for the exact initial 
state ipo(x) = xe _10 ' x ''. 
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FIG. 13: Probability spectrum with ^o(^) = xe 10 \ x \. Parameters t = 0.1, n = 100. 
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FIG. 14: Probability spectrum with i)o{x) = x\x\e 10 \ x \. Parameters t = 0.1, n = 100. 

The exact energy is E = —k 2 /2 = —50, and the agreement is satisfactory. Fig. 14 shows 
the energy spectrum with initial state ipo(x) = x\x\e~ 10 ^, which contains excited states. 

The spectrum shows another bump around E ~ —10, which corresponds to the first 
excited state with energy E x = —k, 2 /8 = —12.5. Fig.15 shows the projected wave function 
corresponding to E = —51.05 of Fig. 13. 
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FIG. 15: Projected eigenfunction of the ground state. Solid line shows the exact wave function, 
(a) Real part, (b) Imaginary part. 
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The phase of the wave function is set to real at the maximum amplitude (\x\ = 3/2 5 ). The 
agreement seems fairly good, although the mixture of the imaginary part is not negligible. 

Fig. 16 shows the average result of 10 random initial states. In this case, only the ground 
state (E ~ —50) can be seen. This is because excited states are accumulated near E ~ 
in the Coulomb potential, and positive energy continuum states might contribute to fill the 
whole energy range due to the periodicity. 

0.2 
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FIG. 16: Average of 10 random initial states. 

IV. SUMMARY 

We have explicitly constructed quantum circuits and carried out simulations of typical 
one-dimensional Schrodinger equations, i.e., harmonic oscillator, square- well and Coulomb 
potential. We have made quantum circuits in such a way that they consist of only single- 
qubit and two-qubit operators and do not require ancillary qubits to calculate the potential 
term. Therefore they are simple and easy for implementation. With eight qubits (4 work 
qubits and 4 simulation qubits), our simulations could obtain reasonable outputs compared 
with the exact results. It is found that exactly symmetric mesh points should be employed 
for the symmetric potential, and the initial states should be prepared deliberately. 
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